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The purpose of this paper is to investigate several analytical methods of solving 
hrst passage (FP) problem for the Rouse model, a simplest model of a polymer 
chain. We show that this problem has to be treated as a multi-dimensional Kramers’ 
problem, which presents rich and unexpected behavior. We hrst perform direct and 
forward-hux sampling (FFS) simulations, and measure the mean hrst-passage time 
t{z) for the free end to reach a certain distance away from the origin. The results 
show that the mean FP time is getting faster if the Rouse chain is represented by 
more beads. Two scaling regimes of t{z) are observed, with transition between 
them varying as a function of chain length. We use these simulations results to 
test two theoretical approaches. One is a well known asymptotic theory valid in the 
limit of zero temperature. We show that this limit corresponds to fully extended 
chain when each chain segment is stretched, which is not particularly realistic. A 
new theory based on the well known Freidlin-Wentzell theory is proposed, where 
dynamics is projected onto the minimal action path. The new theory predicts both 
scaling regimes correctly, but fails to get the correct numerical prefactor in the hrst 
regime. Combining our theory with the FFS simulations lead us to a simple analytical 
expression valid for all extensions and chain lengths. One of the applications of 
polymer FP problem occurs in the context of branched polymer rheology. In this 
paper, we consider the arm-retraction mechanism in the tube model, which maps 
exactly on the model we have solved. The results are compared to the Milner- 
McLeish theory without constraint release, which is found to overestimate FP time 
by a factor of 10 or more. 
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I. INTRODUCTION 


Rouse model is the simplest stochastic model of polymer dynamics, where almost every¬ 
thing can be solved analytically. The polymer chain is modelled with a set of beads with 

coordinates Hi and friction ,^o experiencing Brownian motion and connected into the chain 

3k^T 

by a set of harmonic springs with spring constant k = ——— where ks and T are Boltzmann 

¥ 

constant and temperature and b is the statistical segment length, respectively. Thus, the 
equation of motion of a given bead i is 

^ {R^+l + R^-l - 2R.)+m; im) = o; = 2kBnob^At-^') 

( 1 ) 

where the hrst term on the right hand side represents the connectivity forces and the second 
term is random white noise force. This equation is valid for all beads i apart from the two 
end beads. The equation for the end beads depends on the problem - they can be hxed, 
free, or have a constant force acting on them. In this paper we will consider chains with one 
end hxed at the origin (Rq = 0) and the second end — free. 

Most experimental quantities of interest, such as stress or end-to-end relaxation functions, 
can be calculated analytically for Rouse model. This is done by transforming the system of 
coupled eq|T| into a set of independent equations for the eigenmodes of this system {Xp}, 
called Rouse modes in this context. Equation for each mode can then be solved indepen¬ 
dently, and all physical quantities are expressed as combination of these Rouse modes and 
their correlation functions. 

A notable exception from this is the hrst-passage (FP) problem (often called Kramers’ 
problem), which is still not solved analytically even for the Rouse model. The FP problem 
requires calculation of the average time that the free end of polymer chain reaches a particular 
point for the hrst time after starting from some specihed initial set of conformations (e.g. 
those in equilibrium state). This problem appears in many areas of polymer physics, ranging 
from biophysic^ to chemical reaction^^ to polymer rheologj^^. It is also related to protein 
folding, although in proteins all amino acids must fold into correct places, not just the end 
groups. In this work, we shall concentrate on a particular example from rheology of polymer 
stars, but we believe many of our hndings will be applicable to other scientihc areas as well 
as to general FP problems in multi-dimensional systems. 

In polymer rheology, the dynamics of polymer melts made from long chains can not be 


2 





described by the Rouse model because the chains get entangled with each other, which 
means that their interactions are much more complex than those reflected by eq|g This, 
however, can be approximately solved using the tube model concept: one assumes that each 
chain is moving parallel to the contour of the tube (made by other chains) and the motion 
perpendicular to the tube is suppressed. In order to relax imposed deformation, the chain 
must wiggle out of the tube into the newly created tube which does not carry the memory 
about deformation. If we accept this concept, the problem of stress relaxation in entangled 
polymer melt is reduced to the first-passage problem because the stress associated with 
each tube segment gets relaxed when the chain end reaches it for the first tim^. Since 
the motion of the chain along the tube is assumed to be unaffected by entanglements, one 
usually assumes that the Rouse model eqj^ is valid if we take Ri to be the positions of 
monomers along the tube axis. In other words, in order to solve most problems in entangled 
polymer dynamics, one has to start by solving FP problem for the one-dimensional Rouse 
chain, which is exactly the focus of this paper. We note that entangled polymer dynamics is 
greatly complicated by the fact that the tube itself is not fixed in space due to the motion of 
other chains surrounding the target chain, a phenomenon which is called constraint release. 
We will not address this problem in the current paper, or in other words we will assume that 
surrounding chains are much longer than our target chain and therefore provide permanent 
network of entanglements. 

Before plunging into the depth of the problem, we briefly describe the Milner-McLeish 
theor}^^ currently used in polymer rheology of stars. The theory assumed that the Rouse 
chain inside the tube can be replaced by one bead attached to the origin through a harmonic 
spring. The spring constant was then chosen to be /c = ^^|^such that the average squared 
spring length is equal to the mean-square end-to-end distance of the polymer chain (where N 
is the number of moving beads in the chain). The bead friction was chosen to be ^eff = y'Co, 
i.e. to carry half of the friction of the whole chain. Authors argued that the reason for this 
1/2 factor is that the average bead only needs to travel distance L/2 if the end is to travel 
distance L. The FP problem of the one bead model has an exact solution (Kramers formula) 


TiZ = 


M r 

ksT Jq 


dx exp 


U{x) 


exp ( — 




k^T 


u'{z)y K"(o) 


exp 


ksT 


( 2 ) 


where U{z) = 


SkeTz'^ 

2^2 


The approximation is valid if U{z) S> ksT. This leads to the 
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following expression for the mean FP time 


Tmm{z) 


71 - 5/2 y/Nh / 32;2 \ 


( 3 ) 


where tr = —--is the longest relaxation time of one-end hxed continuous Rouse chain, 

which is 4 times larger than that of the chain with two free ends. 

In this paper we shall hrst perform direct and forward-flux sampling simulations of the 
FP time of one-dimensional Rouse chain and compare them with the Milner-McLeish as¬ 
sumption (section 2). A very signihcant disagreement will motivate us to look at the exact 
asymptotic solution of a multi-dimensional Kramers problem in section 3, which we will hnd 
methodologically useful but invalid for the realistic extensions of Rouse chain. In section 
4 we will develop a new FP theory for the multi-dimensional Gaussian process, which will 
result in correct scaling behaviour but an incorrect prefactor in the intermediate regime. We 
will hnally combine simulation and analytical results into a simple formula for the FP time 
of discrete and continuous Rouse chains for any large extension where U{z) S> ksT. 


II. DIRECT AND FORWARD FLUX SAMPLING SIMULATIONS 

In order to set the scene for further calculations, and to evaluate the Milner-McLeish 
ansatz, eqj^ we hrst perform computer simulations with variable number of beads rep¬ 
resenting the Rouse chain. The bead friction temperature /crT and statistical segment 
length b are set to unity in the simulations without loss of generality, which sets the length(6), 
/ksT) and energy(fcRT) scales. Direct simulations of eqj^were performed using the 
predictor-corrector methocP. Detection of hrst-passage events was improved by computing 
the probability of unobserved events in each time stejP. Even if no crossing is observed 
before and after a time step, there is still a possibility that trajectory has crossed an inter¬ 
face during the step. If qi and q 2 are the distances between the interface and the reaction 
coordinate before and after the time-step, such probability is 

Pcross — exp( 

where D is the short time diffusion coefficient of the reactive coordinate, in our case, of 
the last bead: D = ksT/^Q. A uniform random number w on [0..1] interval was generated 


4 







and the unobserved event registered if w < Pcross- With this algorithm and the predictor- 
corrector integration we were able to use time step At = 0.05, with smaller steps giving 
almost identical results. 

Direct simulation results for the mean FP time are plotted in Figjljwith solid lines for 
different numbers of beads N. The horizontal axis in both plots is the distance from the free 
chain end to the hxed end, 2 ;, normalized by the root-mean-square end-to-end distance that 
s = z/{N^Ab). In a continuously simulation, when the free end last crossing Sq = 0 reaches 
s > 0 for the hrst time, its time cost is recorded as the FP time for s. Figj^a) shows decimal 
logarithm of the mean FP time. Clearly, the time grows very fast with s, approximately 
as exp(3s^/2) as expectecP, and the simulations can only carry out the measurement up to 
r(s) ~ 10^ or so. FigQb) shows the same data in reduced coordinates r(s)sr^^ exp(—3s^/2). 
Here for clarity we divided by the trivial factor exp(U{s)/kBT) = exp(3s^/2) predicted by all 
theories, and the typical fluctuation time tr, and also multiplied by the expected prefactor 
scaling s as suggested by eqi Such a renormalized plot brings all data within one decade 
in vertical axis scale and allows clear comparison between theories and simulations. In 
particular, the Milner-McLeish result, eqj^ is constant in this representation and is shown 
by dot-dashed line. Already from the direct simulations, it is clear that simulation results 
are signihcantly faster than the Milner-McLeish prediction. 

In order to extend simulation results to longer times and facilitate detailed theory ver- 
ihcation and calibration, we also performed forward-flux sampling(FFS) simulation of the 
same model. FF^^^ is a simple method to compute FP times of unlikely events by splitting 
the phase space of the system by n -|-1 interfaces dehned such that the starting points of the 
trajectory are on the left of 0-th interface, and the hnal point — on the last, n-th, interface, 
as shown in Figj^ Besides, the interfaces must be dehned such that every trajectory has 
to pass all interfaces in order to get to the reactive state. In our case, the position of the 
i-th interface can be simply dehned by the position of the free end Rn = A*. We found 
that Ao = 0, Aj = (1 -|- 0.25 * (i — l))N^Pb, i > 1 gives the most accurate results, avoiding 
systematic errors due to very small interface distance and large statistical errors due to large 
interface distance. 

The simulation then proceeds in two stages. In the hrst stage, we run one long simulation 
for time Tq and count the number of crossings, A"o, of the hrst interface Ai by the trajectories 
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FIG. 1. (a) The logarithm of FP time r(s) versus s for FFS simulations (dots) and direct simulations 
(solid lines), (b) t{s)st^^ exp(— 3/2s^) versus s for FFS simulations (circles) and direct simulations 


(solid lines), the dashed lines are the prediction of Equation 31 Milner-McLeish theory is shown 
by the red dash-dot line. 


No 

J-0 


which last crossed interface Aq, rather than Ai. These results are in the attempt frequency 

( 5 ) 

Besides counting the crossings, we store the full chain configurations at the moments of these 
crossings. 

In the second stage, we run many short consecutive simulations for interfaces 1 to n — 1 in 
sequence. For the interface A*, the simulation starts from the stored points on the interface A* 
(selected at random from the database) and finish when they either reach the next interface 
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FIG. 2. The sketch of the direct and FFS simulation. The Rouse chain with one end fixed do one¬ 
dimensional Brownian motion, z can be normalized into s for the comparison between different 
chain length. The interface definition is only for FFS simulation. 

A(successful run), or go back to the 0-th interface (unsuccessful run). The fraction of 
successful runs Ni/Mi gives an estimate of the probability to progress from one layer to the 
next, P(Aj+i|Aj), where Mj is total number of runs from layer i, and iVj is the number of 
successful runs. Thus, the mean FP time is given by 

r(A„) = --. (6) 

The value of Mi has a decisive effect on the statistical error of the dual outcome, with 
the best strategy to increase Mi for higher energy barriers between the layers to ensure an 
approximately constant iVj. A simple way to determine Mj is to run a few simulations with 
smaller Mi and get the rough ratio of P(Aj+i|Aj), and estimate Mi for an expected A,. Ref.^^ 
recommends selecting interface distances such that P(Aj_|_i|Aj) > 0.3. Our selection satishes 
this criteria. By running a quick simulation for N = 1, the proper Mi can be obtained. 
Using the same Mi and the same distance dehned by s for A = 1, a proper ratio P(Aj+i|Ai) 
for larger A is also guaranteed since the P(Aj+i|Ai) increases with larger A. 

The mean FP time r(s) for various normalized extension length s are presented for 
different chain lengths A in Fig{^ The FFS result is the harmonic average of FP times from 
independent runs. More discussion on the averaging method and the error can be found in 
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Appendix A. Comparing with direct simulations, two methods disagree slightly for s < 1.5. 
In this region, the FP time given by FFS is inaccurate since the energy barrier is below 
S.bksT, different FP events are not independent from each other and therefore distribution 
of FP times is not single exponential. In the region of s > 1.5, two simulation methods are 
consistent with each other. FFS method is able to predict the FP time till s = 5.5, and the 
chain length up to iV = 128. In the normalized plot in FigQb), all curves develop from 
a negative slope to a plateau after exceeding some transition length St- The slopes of the 
curves from the peak to s* increases with increasing N. In the mean time, the transition 
length St also increases. One hnds that the result differs from the Milner-McLeish theory 
signihcantly. When increasing N, the FP time becomes much shorter then their prediction, 
leading to the difference of a factor of 10 at s = 3 and N = 128, and even bigger for larger 
s and N. This shows conclusively that the one mode assumption of the Milner-McLeish 
theory is inadequate and better theory must be developed. Note that this discrepancy is 
much bigger than the 20% reported by Vega et.al.ll^. It suggests that the multi-dimensional 
character of the end bead evolution is essential, motivating the next two sections. 


III. EXACT ASYMPTOTIC THEORY 


In this section we step away from the one-mode approximation of the Milner-McLeish 
theory and apply exact asymptotic Kramers theor}ff2HISl|;o the full problem in A-dimensional 
Rouse modes space. We hrst briefly recall the theory here. Consider a Brownian particle 
with friction ^ diffusing in a multi-dimensional space of dimension d, subject to the potential 
U{R). In the over-damped regime, its probability density distribution obeys the 

Smoluchowski equation 

= v(^(R,t)VV(R) + kBTVij{R,t)) (7) 


which can be divided into two equations 


d'lp 

J{R,t) 


-VJ 


where J{R, t) is the current at point R. The hrst equation is called continuity equation and 
the second — Pick’s law. 




Since the equilibrium solution of such equation is given by Boltzmann distribution 
ipeqiR) ~ exp {—U{R)/kBT), it is convenient to change the unknown function to 


0(i^, t) = t) exp{U{R)/kBT) 


such that (j){R, t) becomes constant in equilibrium. Then eq@ becomes 


{ f U(R)\ 3<l>(R,t) 

kBT \ kBT ) dt 


V ( exp 


( uiR)\ 
V kBT J 




and the expression for the current 

= -^^exp V0(i^,^) 


( 8 ) 


The Kramers method of hnding the mean FP time of escape is to consider a system where 
particles are injected at the origin and deleted from the system when they reach absorbing 
surface and hnd a steady state solution of such system. In this case the average escape 
time will be given by the total number of particles divided by the rate of injection 


T = 


1 

Jin 


tlj{R)dR = 


1 

Jin 


((>{R) exp 


( uiR)\ 

V kBT J 


dR 


(9) 


where integration is over the whole space available to particles and ipiR) and 4>{R) without 
time dependence are steady state values. 

An asymptotic result as T —)■ 0 was obtained and proved rigorously by Meerko\fl^^, 
based on earlier papers of Kifei®. Here we show a simple intuitive derivation which is not 
available in the original papers. 

The crucial assumptions are that as T —>■ 0, (j){R) becomes constant everywhere apart 
from a thin layer near the absorbing interface, and in this thin layer the current is perpen¬ 
dicular to the surface. Let’s call x the direction orthogonal to the interface and Y = {j/j} 
— all other directions. Then in the layer near absorbing surface we can write 


d<p _ ^UY) fUiR)\ d<P _ 
dx kBT """"PVfcBTj’ % 

where Jx{Y) is the current near the absorbing surface. The solution can be written as 
(l){x,Y) = [ exp{U{x',Y)/kBT)dx' ^ exp{U{xs,Y)/kBT) 


( 10 ) 


( 11 ) 
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where Xs is the value of x and the surface. Because we want 0(a;, Y) to be constant far away 
from the surface, we require that 

<P{x, Y) = exp(f/(x„ Y)/kBT) =C ^ UY) = exp(-[/(a:„ Y)/kBT) 

U'Axs^Y) i 

The last equation is sufficient to calculate the total current through the surface, which must 

be equal to the injection rate 


Jtot — Jin 


.h{xs,Y)dY 


C 

1 


t/'(x„l-)exp 


U{xs,Y] 

kBT 


dY 


( 12 ) 


Let E to be the Hessian matrix of the potential U{x, Y) under the new rotated coordinates 
and we partition E as follows. 


E 




(13) 


where E^x is the entry of E at the hrst row and the hrst column corresponding to x-direction, 
A is a column vector with iV —1 elements and E' is (iV—1) x (iV —1) square matrix. Since the 


exponential term in eq. 12 can be treated as the distribution of multi-variables Y conditional 
on a: = Xg, this conditional distribution is multivariate normal {Y\x = Xg) ~ A/’(0, S) where 
covariance matrix 


i: = E'- AE-^A^ 


(14) 


such that 


det(^^) = det(ii^a; 3 ;) det(S) 


(15) 


Suppose the minimum of the potential on the absorbing surface is reached at {xg,Ys). 


The integral of eq 12 is dominated by the potential in the area around {xg,Yg) and can be 
approximated as 


Ji 


CU'x{xg,Yg) 


eM-U{xg,Yg)/kBT) 


(27rA:BT)(‘^-i)/2 


(16) 


The mean FP time is then the ratio of total number of particles Ntot to the total current. 
The major contribution to the number of particles comes from the bottom of the potential, 
where ((>(i?) is approximately equal to constant C\ 
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= U{R)exp{-U{R)/kBT)dR 


{2TikBTY/‘^ 


C 


(17) 




Vdet(Ai 

where Ajj is the Hessian of the potential computed at the minimum of the potential in the 
original coordinates. 


Combining eqs, 16 and 17, we arrive to the final expression, similar to the one from 
Meerkov’s paper 


N. 


Hot 

^ - - OV 




1 27r/cB7"det(S) 


exp{U{xs,Ys)/kBT) 


(18) 


Jin ^^s) V det(Ajj) 

We note many similarities with the one-dimensional result. In the multi-dimensional version, 
the first derivative of the potential at the arrival point Xg is replaced by the hrst derivative 
along the normal to the surface. The second derivative at the bottom is replaced by the 
determinant of the Hessian, both being proportional to the volume available to the particles 
at the bottom of the potential, or the volume of the area where U{R) < kBT. The new 
term det(S) has a meaning of channel width available to the particles as they arrive to the 
most probable absorbing point {xs,Ys). 


Now we will apply the result of eqj^ to 1-D discrete Rouse chain. If one end of the 
Rouse chain is attached to Rq = 0, the transformation to eigenmodes Xp(or Rouse modes) 
and back is given by 

^ D 1 / 2 ). D _ o 1 / 2 ) 

y Rj sin j\f 7-^2 2 Xp Sin 7^?- , /r> 

p=l 


N + 1/2 


2=1 


iV + 1/2 ’ * ^ ^ N + 1/2 

p=i 

Since one end is hxed at origin, the end-to-end vector of the Rouse chain is just 


N 


Rn ^ ^ cipdfp, 

p=i 


( 20 ) 


where 


ap = 2 sin 


7iN{p — 1/2) 

N + 1/2 ' 


( 21 ) 


The equations of motion for Rouse modes are then. 


dXp 

dt 


=-f3pXp +fp{t); p=l,2,---,N 


2kBT 


(fpit)) = 0 , iUYfYt')) = -^5{t - 

sp 


Jpq 


( 22 ) 

(23) 


11 

















where fp{t) are the random white noise and f3p are dehned as the ratio between the spring 
constant Kp and the friction coefficient ^p which are 


Kr 


2AkBT{N + 1/2) . 2 7r(p-l/2) 
siri 


62 


2{N + l/2) 


^p =2{N + l/2)eo 


(24) 

(25) 


such that 


I3p 


2AkBT{N + 1 / 2 ) ^.^2 - 1 / 2 ) 




2{N + l/2) 


(26) 


As shown in eq,25, the values of ^p are the same for all p, therefore we denote them as 


= ^p,P = ' ,N in the rest of the paper. Here, we also introduce two useful sums used 

frequently for the Rouse chain in the rest of paper, 

N 

^aJ=2]V + l (27) 


p=l 
N 

p=l ^ 


al 2N{N + l/2)eo6^ 


(28) 


Because det(Ajj) is independent of the rotation of the coordinates and the potential is 
harmonic in each dimension U{X) = Y/p where ^xPp = is the spring constant 

of the potential in the p-th dimension, so det(Ajj) = Because the absorbing 

boundary has the expression Yp (^pXp = z, the direction of x is parallel to the unit vector 


{%} = 




v' ! 
V p' 


Since det(Aij) = det(^^) at every point in the coordinate space and 


equals to the product of the fluctuation amplitude in ^(-direction {EY and det(S) according 


to eq 


15 


then \Er 


, det(Ajj) , . , . , . 1 , . r 1 

IS equal to —— 7 =/-, which is nothing but the inverse of the spring 


_ det(S_) 

constant of the potential in the ^(-direction so that 


det(A. 




Yv « 


2 

p —p 


det(S) Y/ 


ql 


p ix+ 


V 2z. 

Xip Pp 


(29) 


The derivative of potential in ^(-direction at is ■ Qp = Yp^xl3pX^ 




Since is the location of the potential minimum on the absorbing boundary 




, we have 


t/i(x„y,)=5^CYA 


OLr, 


zar 


Yp' ^p' /5p Yp' B, 
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( 30 ) 




Pp 


Substituting eqsj^ and into eq|T^ we get 

^0 


t{z) 


ksT 3v/6 ^ 

2 


-^ 3 / 2^3 fUiXs.Y 
exp 


327r 


sm 


TT 


4(Ar + l/2)/ ^ 


ksT 

N^/% 


-tr exp 


where the Rouse time of the discrete chain is = ^ = 


- J_ - gpfe" 


ksT 

— 2 TT 


sm 


4(W+l/2)' 


(31) 


These asymp¬ 


totic results for different N values are shown in Figj^b) by dashed lines together with the 
ones obtained from FFS simulations. In the limit of large N, the result in eq,^ can be 
simplihed as 


r(z = 


7 r 5/2 Tr y/Nb ( 3^2 ^ 




(32) 


2^6 ^ 

which is ^ times smaller than the Milner-McLeish result, eqj^ Interestingly, this factor is 
exactly the ratio of friction of one bead to the effective friction used by Milner and McLeish, 
which means that the asymptotic result can be obtained by solving one-dimensional problem 
with spring constant k = 3kBT/{Nb‘^) and the friction of one bead. This means that in the 
limit of zero temperature or very large extensions the friction of other beads must be 
neglected. 

One of the most important conclusion is that this result can not be expressed as a function 


of reduced variable s and t/tr only because it has an additional factor of A; 


r(s) 1 7 r ®/2 X 


exp 


2 


(33) 


Tr N 2\/6 s 

Thus, it becomes arbitrary small in the limit of continuous Rouse chain. This does not seem 
to be physical, and a more detailed theory of the next section will reveal the reasons behind 
such behavior. 


IV. PROJECTION ONTO MINIMAL ACTION PATH 

In mathematics, a widely known approach to FP problem is the Freidlin-Wentzell 
theorjff^, which operates with the probability of the Brownian paths in terms of an ac¬ 
tion functional. It uses the fact that the most probable path is given by minimizing the 
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action functional associated with the system. Sophisticated numerical techniques have been 
developed to hnd such transition path and transition rates in complex systemd^S^^^^. In 
this section, we will apply this approach to our system. For simplicity, we hrst consider 1-D 
case. 


A. Minimal action path 
1. 1-dimensional case 


Suppose a particle is injected at x = 0 and t = 0 into a 1-D system with harmonic 
potential U{x) = where the spring constant k = (3^, ^ is the friction coefficient of the 

particle. This is one of the simplest stochastic process called Ornstein-Uhlenbeck process. 
The particle fluctuates in the potential and is deleted from the system at x = z and t = t*. 
In order to identify the trajectory of the particle, a small time interval for recording the 
position of the particle is set to At such that the trajectory consists of n segments and 
corresponding n -|- 1 coordinates where n = t*/At. In other words, the trajectory of the 
particle can be somehow treated as a “bead-spring chain” with n + 1 beads where two ends 
hxed at X = 0 and x = z, respectively. 

Since Ornstein-Uhlenbeck process is Markovian, the chain segments are independent from 
each other. The probability of the chain having conhguration {r} is simply the product of 
the transition probability density associated with P(rj+i|rj), 

P{{ri,r2, ■ ■ ■ ,rn+i}) = P(r„+i|r„)P(r„|r„_i) ■ ■ ■ P(r2|ri)P(ri) 


The transition probability for Ornstein-Uhlenbeck process is well known 


P(rj+i|ri) = 

where D = ksT/^. Then 


1 


v/27rP(l - e-2/3At)//5 


P({r}) = 




27rP(l - e-2^^0 


n/2 


exp 


exp ( — 




2D{1 - e-2/3At) ^ 


/9(ri+i-rie 
2D{1 - 


^(ri+i - TiC 


(34) 


Once the number of segments n and t* are chosen. At is constant such that the prefactor in 
front of the exponent term is just constant which only contributes to the normalization of 
the probability. Thus, 

r d ” 1 

P({r}) ~ exp ^ - 2/^(1 _ e- 2 / 3 At) ~ (if At is small) 
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' exp 


= exp 


1 ” Ar 2 


AD 


2 = 1 


At 


+ 2riAri/3 + r‘f/3‘^At) {let At —)■ 0) 


~4D I ^ r{t)flfdt 


( 35 ) 


diT 

where r{t) = 

dt 

According to eq,35, finding the most probable trajectory {r} is the same as finding the 
trajectory which minimizes action A({r}) = J [r{t) If we define the Lagrangian 

of the system as the integrand of the action L{r, f, t) = [r{t) + r{t)/3Y, the action will have 
an extremal if the Euler-Lagrange equation is satisfied along the trajectory 


d dL 
dt dr 


dL 

dr 


(36) 


Since the second derivative of the Lagrangian over space and momentum are both positive, 


the action has a minimum along the trajectory. Eq,36 leads to an ordinary differential 
equation for the minimum action trajectory 


with the general solutions 


r{t) = /3‘^r{t) 


r{t) = Ciexp(—/Sf) + C2exp(/3f) 


(37) 


(38) 


where Ci and C 2 are arbitrary constants. 

In order to capture the correct optimal path, the selection of a large enough time interval 
[0, t*] is essentiapZES. If the absorbing boundary is far from the origin so that the energy 
differences between points on the boundary and origin are much larger than ksT, the optimal 
path is not sensitive to the value of t* as long as t* is reasonably large. To simplify the 
problem further using the reversibility of equilibrium dynamics, we can assume that the 
particles are injected into the system at the absorbing boundary x = z a.t t = 0 and will 
eventually arrive in origin such that r(0) = ^ and r(+cx)) = 0. Because the second term in 


eq,38 is divergent, r{t) will drop to origin if and only if C 2 = 0. Therefore, the solution with 
the specific boundary conditions is 


M{t) = r{t) = zexp{-/3t) 


(39) 


which is termed as minimal action path(MAP). Note that the parameter t here is not the 
real time but a convenient parametrization. 
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2. N-dimensional case 


In A^-dimensional case, suppose the absorbing boundary is defined by the expression 
J2p^pXp = where p E {1, 2, ■ ■ ■ , N} and Op are real numbers. To obtain MAP M{t) we 
first notice that the action in multi-dimensional case is just a sum of actions in individual 


dimensions, and therefore eq,37 and its solution eq,39 apply 


M.p(t) = C„exp(-/?„i) 


(40) 


To find the unknown coefficients Cp, we use two initial conditions. First, the starting point 
must lie on the absorbing plane 

N 

apCp = z (41) 

p=i 

And second, near the absorbing plane the direction of MAP must be perpendicular to to 
absorbing plane. This is because on small length scales the dynamics is dominated by the 
thermal fluctuations and therefore potential can be neglected, which means that the particle 
will find the shortest path to the absorbing plane. Combining these two conditions leads to 
the solution 


Mp{t) = 


ZOr. 


exp(-/3pt) 


(42) 


l^p Y.p' 

It is easy to check that the arrival point M„(0) = {zap/(3p)/ corresponds to the 

minimum of the potential on the absorbing surface. 

As an example, Figj^ shows the MAP between the origin and the location of minimal 
energy on the hyperplane x -f y = 8 in the system with potential U{x,y) = /2 + lOj/^/2, in 

which the curved lines are the contour plots of the potential. Note however that for smaller 
z the actual minimum action path will deviate from our predictions since the action derived 


from eqj^ was computed without the effects of absorbing boundary. 

The key idea of our method is to map the A^-dimensional Kramers process onto a 1-d 
problem along the MAP. The effective potential along the MAP is defined as a function of 
the distance I away from the origin along the path such that the mean FP time for particles 
reaching the absorbing boundary can be derived from Kramers’ solution eq§ 

Since we are interested in the FP time when the particle hits the absorbing boundary, 
the trajectories should not go across the boundary. If the absorbing boundary is far away 
from the origin, the portion of the trajectories exceed the boundary is negligible. However, 
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FIG. 3. MAP(eqj42[) between the origin and the location of minimal energy on the boundary in 
the system with potential U{x, y) = x^/2 + lOy^/2. 

if it is close to the origin, the particles would hnd shorter path reaching the boundary so 


that the MAP of eq,^ is not a good approximation of the most probable trajectory in the 
real system. 

Without loss of generality, we assume /3j < f5j ii i < j so that the slowest mode is Xi. In 
Kramers’ solution, the inner integral would be dominated by the minimum of the potential 
for large values. As shown in Figj^ the MAP follows the x-axis near the origin which is 
the slowest mode such that the potential on the MAP near the origin is nothing but 
where is the smallest spring constant among all Thus, using eqig 

I ^nksT 

Mx Jo 


t{z) 



^l{z) 


kuT 




(43) 


Because the MAP is given by eq,42, the distance l{t) between the origin and M{t) along 
the MAP can be obtained by 


and 


m = 


di 

dt 




dt J 





(44) 


(45) 


( 46 ) 
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Clearly, MAP only gives the information about the most probable trajectory of the par¬ 
ticles. Along the MAP, the particles fluctuate in a relatively narrow channel with a certain 
width which depends on the potential landscape and the direction of the MAP. Since the 
MAP has a sharp turn approaching the absorbing boundary, it is reasonable to believe 
that the width of the channel would also change when approaching the boundary. We will 
investigate these effects in the next section. 


B. Conditional entropy 


In this section, we will investigate how does the particle fluctuate perpendicular to the 
MAP. It is reasonable to believe that the distribution of particle positions in the direction 
perpendicular to the MAP is Gaussian if the particle is far away from the absorbing bound¬ 
ary. This assumption provides a simple approach for obtaining the entropy S{t) of transverse 
fluctuation along the trajectory. However, the distribution might deviate from Gaussian as 
the particle approaches the absorbing boundary. To the hrst approximation, the density 
distribution of the particle is proportional to the equilibrium distribution in the absence of 
the absorbing boundary, so that the entropy 


S{t) = —ks / P{x\x G H(t))\nP{x\x G H{t))dx 


(47) 


- E -^ 


), H{t) = 6{q{t) ■ {x - 


where P(x) ^ Peoix) = -7777- ^ ^ ■ / „ n, ^ 

M(t))) is the hyperplane perpendicular to the MAP and passing through M(t), q(t) = 

tangent vector along the MAP and is the distance from the 
origin(t = -|-cx)) along the MAP. 

Thus, 


P(x\x G H{t)) 


Peq{x)5{q{t)-{x-M{t))) 
f Peg(x)d(q(t) ■ (x - M(t)))d 


X 


(48) 


Because the Dirac ^-function can be represented by <^(A) = ^ the denominator 

can be easily calculated as ("27r/c^T ^ ) exp ( -• Substituting 


m eq 


48 


eqs,24, 25, 26, 48 into eqj47] gives 
S{t) =kB { In 




ksT 


2kBTxe 


\n 2kBT7re^ 


f^pix 


(49) 
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We can see that in the second term in eq 

PpiX ^ 


49 


is nothing bnt the amplitnde of particle 


flnctnation in the direction q{t). At t = +C)0, qp{t) are 0 for all p except p = 1 since the 
MAP follows the slowest mode near the origin, as illnstrated in Figj^ so that 


S{+oo) = ks In , 

p>2 


2kBT7ie 

(^pix 


(50) 


which is exactly the entropy of {N — l)-dimensional mnltivariate normal distribntion. We 
are only interested in the entropy difference AS{t) between points on the MAP and the 
origin, which is given by a very simple expression, 

AS(t) = S(t) - S(+oo) = ^ In 5^ (51) 

p 


C. Effective potential along MAP 


The effective potential along the MAP can be dehned as Ueff{t) = U{M{t)) — TAS(t). 


Snbstitnting it into eqj^ we get onr central result for the mean hrst-passage time, 

r+'Ao o 2 , 

r(.P=V^l ^y-F(lV.l)exp(— 

where 


(52) 


F(N,t) 


N(xb^ 



^ exp{-2/3pt) 
Pp 


(53) 


We note that F{N, t) is related to the simple correlation function (p{t) = {R{t)R{0))/{R‘^{t)) 
with R(t) = Yhp^p^pi'^) without absorbing boundary by F{N,t) = (p{2t). In addition, we 
consider the entropy corrections at every point on the MAP instead of only at the absorbing 
boundary in section.3 and some other literatured^, which, as far as we known, is the hrst 
try. 


D. Discrete Rouse chain 


Now we will apply our theoretical prediction to the discrete Rouse chain and t{z) in eq.52 


is calculated using the parameters in eqs,21 and 26 After changing variable t' = t/rn, the 
mean FP time for the chain end is 


t{z) = 


f*+oo 


VNb 


\/^WOexp(^F(iV,F))cit' 


(54) 
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where 


F{N,t') 


N 2 / 7 rAf(p-l/ 2 ) 

V ^ ^ exD f - 2( 

„• _2/7rfp-l/2) N ^ 


2]V(iV + 1/2) ^ 


Sin 


2 (jV+l/ 2 ) Y 


‘0 


p=i V 2 (Ar+l/ 2 )>' — 4 (/V+l/ 2 ) 

If ^ 1, F{N,t') has three different regimes as follows(t' is in the unit of r^j), 


(55) 


1 - 

t' < 1/N^ 

1 _ /32t' 

y TT^ ’ 

1/N^ < F < 1 

8 -2t' 

, 

t' > 1 


where the last regime is only contributed by the slowest mode and details of the derivation 
can be found in Ref.® 

For the reduced distance s = used in the previous parts of the paper, the result 
further simplihes to 

/ +00 _ \ 

\/F{N, t') exp y—^F{N, t')j dt' (56) 

Figj^ shows normalized mean FP time t(s), similar to the ones in Figj^b), for different 
number of Rouse beads N. Clearly, we observe two different scaling regimes and the results 
at large s values coincide with the prediction of the asymptotic theory in section 3. With 
increasing N, the range of the first regime increases. 



FIG. 4. Normalized FP time of the finite Rouse chain with different N predicted by MAP pro¬ 
jection theory (dashed lines from top to bottom correspond to N =2,4,8,16,32,64,256,1024). The 
asymptotic theory results from section 3 are shown by solid lines. 


The integral in eq,56 for large s is dominated by small t values such that 

f+o° /3s^ \ 

t{s) ^V^trS\^F{N^ / exp y—[F{N,0) + F\N,0)t)jdt 
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Ci(iV) 

Ns 


r^exp 



C'i(iV) 



71 


4{N + l/2) 


(57) 


We can see that eqj^ is the same as eqj^ so that our theory gives the correct asymptotic 
results of r(s) at large s values in the systems with finite Rouse modes. Figshows Ci(A^) as 


. 5/2 


a function of number of Rouse modes N. If N is large, Ci{N) goes to a constant 


3.57. 



FIG. 5. Ci{N) as a function of number of Rouse modes N. 


In the hrst regime, the integral in eq,56 is dominated by F{N, t) in the range 1/A^^ <t< 1. 
Thus, r(s) can be approximated as following. 


s) Pexp (®d(l - 




tr 


36s3 


exp(3s^/2) 


(58) 
in the 


which is illustrated by black dot-dashed line in Figj^ showing the scaling of 
prefactor. 

After comparing eqj^and eq.^, we see that the transition point of two scaling regimes 
is around s* ~ y/N. In other words, the transition happens at z* ~ Nb where b is the 
statistical segment length. Thus, the prefactor of r(s) will have s~^ scaling until the chain 
reaches its full extension and then it will go back to scaling behavior. 

A simple scaling argument can be used to explain such behavior. The Rouse chain under 
stretching can be treated as a sequence of Pincus blob in the stretching direction, such 
that the structure of the chain inside the Pincus blob is not significantly perturbed by the 


stretching force. The asymptotic theory of section 3, eq,33, predicts that the FP time should 
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scale as tr/(N s) exp ^3/2s^^, but this formula is only valid in the limit of fully stretched 
chains, s > \fN. If we chose the level of coarsegraining to be the Pincus blob, we will satisfy 
this condition. In other words, to be physically meaningful, equation should have the 
number of Pincus blobs instead of N in the prefactor. Thus, at particular z, the number 
of beads one must represent the Rouse chain with should be equal to the number of Pincus 
stretching blobs. This number is easily computed from the condition that the sizes of all 


blobs, each containing g monomers, must add up to ~ giving the number 

of blobs 


AT _ ^ _ 
^blobs 


Nh^ 


= s 


Using Nbiobs in eq 33 instead of N changes the scaling of r(s) to s 


-3 


(59) 

This argument also 


explains the transition to asymptotic regime r(s) ~ s which occurs when the number of 
Pincus blobs reaches N, giving s* ~ a/N as discussed above. 


E. Comparison with FFS 

The comparison of r(s) between Forward-Flux sampling simulations (symbols) and our 
theoretical predictions (lines) is shown in Figj^ The number of Rouse beads N used in the 
simulations are 1, 2,4, 8,16, 32, 64 and 128 from top to the bottom. The agreement between 
results with iV < 4 is reasonably good. With iV = 1, our theory gives the same result as 
conventional 1-D Kramers’ solution so that the prediction from our theory lays on top of 
the simulation results exactly. However, with increasing number of beads, the agreement 
is getting worse in the regime we are interested in. In such regime, our theory gives a 
correct scaling with s, but with the incorrect prefactor which probably due to the fact that 
the process is not Markovian anymore after we project the trajectories onto the MAP. In 
the following section, we will provide an empirical expression to overcome the defect in the 
theory. 


F. Combining analytical and simulation results 


Since we know both the exact asymptotic solution of r(s) at large s values and the one 
in the intermediate regime, an empirical expression 

, , fC,(N) , C2(JV)^, 
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FIG. 6. Comparison of r(s)s/ exp(1.5s^)/rK between direct FFS simulations (circles) and numerical 
calculation based on eq|52|(lines) with different number of beads(A^ = 1,2,4,8,16,32,64 and 128 
from top to bottom). 


can be proposed to recover the results in both regimes and to approximate the integral in 


eq,56, where C'i(iV) is given by eq 57 




FIG. 7. (a) Empirical fitting (lines) of theoretical prediction of r(s) (symbols) with different N 
values(A^ = 2,4, 8,16, 32, 64, 256 and 1024 from top to bottom), (b) Fitting parameter C 2 {N) as a 
function of N. 


Fig0 a) shows the fitting of our theoretical calculations by using eqj^ where C 2 {N) is 
a htting parameter, and Gi(iV) is hxed by eq. 57 The htting parameter C 2 {N) shown in 
Figj^b) increases with increasing N and eventually reaches a plateau. The relative fitting 
error, defined as the maximum of - approx each curve, is smaller than 5%. 

exact 

Figg a) shows the htting of FFS results using eqj^ where C 2 {N) is again the htting 
parameter. The ht is extremely good for all chain lengths and for s > 1.5. For smaller s 
neither Kramers’ theory, nor FFS method are applicable. Figgb) shows that the htting 


23 



















FIG. 8. (a) Empirical fitting(lines) of r(s) from FFS(symbols) with different N values(A^ = 

1,2,4,8,16,32,64 and 128 from top to bottom), (b) Fitting parameter C 2 {N) as a function of 

N. 


parameter 6*2 (A^) reaches a plateau at smaller N values comparing to the one in the theo¬ 
retical calculations. The value of 6*2 (A^) at large N is about 3 times smaller than the theory 
predictions. Thus, a combination of theory and FFS simulation leads to a simple analytical 
expression for the FP time valid for N > 10, 


'r(s) = tr 


3.57 1.19\ 



(61) 


V. SUMMARY AND CONCLUSIONS 

In this work, we have studied the first-passage problem for the chain end of discrete 
and continuous 1-D Rouse chains, as a proxy for dynamics of arm-retraction of isolated 
star polymers in a network. In the widely known Milner-McLeish theory, star arms are 
represented by Rouse chains inside their conhning tubes and further replaced by one bead 
attached to the branch point by a harmonic spring. The mean disengagement time of a 
tube segment t{z) is nothing but the hrst-passage time of arm-end reaching it, which can 
be calculated through Kramers’ solution T{z)r^ z-^eMU{z)/kBT). 

In order to check the validity of the Milner-McLeish theory, we carried out direct simula¬ 
tions to collect the first-passage time of Rouse chain extension and found that the mean FP 
time drops signihcantly if the arm is represented by Rouse chain with more beads instead 
of a single bead. Since the large deviations of the Rouse chain happen rarely, an efficient 
simulation method called forward-flux sampling was applied to obtain results for large ex¬ 
tensions. The results from FFS simulations show that the prefactor of mean FP time t{z) 
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has z~^ scaling only for very large extensions, bnt different scaling behavionr t(z) ~ z~^ 
in the intermediate regime. We argue that the large deviations of Rouse chain should be 
considered as a multi-dimensional Kramers’ problem. 

In order to solve the multi-dimensional Kramers’ problem, we hrst apply an asymptotic 
theory, valid in the limit of zero temperature. The solution has many similarities with the 
one-dimensional case. We found that the theory becomes valid only for very large extensions 
~ Nb, corresponding to a fully extended chain, and fails in the physically relevant regime 
VJVb <z<Nb. 


To describe this regime of moderate extensions, a new theory was proposed, which in¬ 
volves projecting the dynamics onto the most probable trajectory between the origin and 
the absorbing boundary, termed “minimal action path” (MAP). We found that in case of 
Gaussian processes, MAP has a simple analytical expression. In addition, we calculated the 
fluctuation perpendicular to the MAP at every point on the path and added such entropic 
correction to the effective potential along the MAP. Then, the mean FP time can be derived 
based on the conventional Kramers’ solution along the MAP. Our theory predicts that the 
prefactor of mean FP time has z~^ scaling before the discrete Rouse chain reaches its fully 
extension and z~^ otherwise. For the continuous Rouse chain, the scaling is always z~^ since 
the chain is never fully extended. Despite of more detailed treatment and correct scaling 
predictions, the constant prefactor in our theory for the hrst regime is about 3 times larger 
than the one obtained from the FFS simulations. This is probably caused by the projection 
of the A^-dimensional system onto a one-dimensional reaction coordinate, which contains a 
hidden Markovian approximation. This hypothesis is supported by the fact that our solution 
for continuous chain, eq.^, is identical to the result of ref.^, although it used very different 
method and the integral expression in ref.^ is different from our result. Since the method 
used in ref.^ was later shown to be suffering from a hidden Markovian approximation!23!^ we 
suspect that our new method did not manage to avoid similar approximation. 


To rectify the disagreement, we combined simulation and analytical results into a simple 
formula for the mean FP time of discrete and continuous Rouse chains for large extension. 
In the future, the analysis of mean FP time of star polymers based on real systems will be 
carried out and the theory will be examinecP^. 


25 



Appendix A: Statistical errors in the Forward-Flux sampling simulations 

Because FFS method works in a propagative manner, the simulation on each interface 
cannot be carried out simultaneously, which restricts parallelization. Meanwhile, it is ex¬ 
pensive to explore the phase space especially in high-dimensional cases. For example, for 
the case of chain length N = 128, the time that the chain evolves from one conformation 
to another uncorrelated conformation is ^ 2230. The two conformations represent only 
two separate regions in a broad phase space, which is far from being enough to provide a 
good distribution on the first surface Ai. In order to improve the accuracy with reasonable 
computational cost, one usually conducts a series of FFS simulations independently from 
different starting conformations, and calculate the FP time by averaging. This idea is equiv¬ 
alent to exploring the phase space from scattered points, although single FFS simulation 
can only explore locally, more simulations will be able to sample the whole space. 



s 

FIG. 9. A comparison between arithmetic mean and harmonic mean for averaging independent 
FFS runs. The simulations are carried out under two conditions, in condition A, expected Ni is 
1000, and there are 100 independent FFS runs, in condition B, expected Ni is 10000, and there 
are 10 independent FFS runs. The inset is the variance for separate conditions. 

For discussion of the systematic and random error, we will consider the case of A = 32 
as an example. The independent simulations are carried on different CPUs, each CPU 
runs equal amount of FFS simulations. For all independent runs, the number of runs from 
each interface, Mj, must be identical, so that their statistical weights are equal. The time- 
step of At = 0.01 is chosen. Fig|^ compares the results from two sets of parameters. Mj 
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is estimated from an expected Ni, which are respectively 1000 and 10, 000. For the hrst 
set, there are = 100 independent rnns, while for the second, there are 10 rnns, snch 
that the time-cost for both sets are eqnal. M, in hrst set is 10 times smaller than in the 
second case, thus the variance is much bigger. By arithmetic mean, the average FP time is 
(r(n)) = T{n), however, a systematic error can be detected after = 4.5, 

where r suddenly goes up and diverges quickly from the average results of the other set, 
indicating that the direct average method is wrong for large ; 2 . In contrast, harmonic mean 


(r(n)) = 


iV. 




(Al) 


provides accurate results for all z. In Figure the FP time using this average method 
converge for both sets. This method allows us to simulate plenty of independent FFS runs, 
whose results can be averaged and are comparable to a simulation with larger Mj, while the 
time-cost in each processor would be much less. 


ACKNOWLEDGMENTS 

This work was supported by the Engineering and Physical Sciences Research Council 
(EPSRC), grant EP/K017683. 

REFERENCES 

^C. Jeppesen, J. Wong, T. Kuhl, J. Israelachvili, N. Mullah, S. Zalipsky, and C. Marques, 
Science 293 , 465 (2001), 

^P. G. de Gennes, J. Ghem. Phys. 76, 3316 (1982), 

^G. H. Fredrickson and L. Leibler, Macromolecules 29, 2674 (1996), 

^A. Kolb, G. Marques, and G. Fredrickson, Macromol. Theory. Simul. 6, 169 (1997), 

®S. T. Milner and T. G. B. McLeish, Macromolecules 30 , 2159 (1997). 

®S. T. Milner and T. G. B. McLeish, Phys. Rev. Lett. 81 , 725 (1998). 

^M. Doi and S. Edwards, The Theory of Polymer Dynamics (Oxford University Press, 
1988). 

®M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Glarendon Press, New 
York, NY, USA, 1989). 


27 




C. Ottinger, Stochastic Processes in Polymeric Fluids: Tools and Examples for Devel¬ 
oping Simulation Algorithms (Springer, 1995). 

J. Allen, P. B. Warren, and P. R. ten Wolde, Phys. Rev. Lett. 94, 018104 (2005), 
Kratzer, A. Arnold, and R. Allen, J. Cliem. Phys 138 (2013). 

A. Vega, J. M. Sebastian, W. B. Russel, and R. A. Register, Macromolecules 35, 169 

( 2002 ). 

I. Kifer, Theory of Probability and its Applications, Vol. 19 (1974). 

^^S. Meerkov and T. Runolfsson, IEEE Trans. Autom. Control 33, 323 (1988). 

^®P. Hanggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990). 

Freidlin and A. Wentzell, Random Perturbations of Dynamical Systems (Springer, 
1998). 

V. Ionova and E. A. Carter, J. Chem. Phys. 98, 6377 (1993). 

^®R. Olender and R. Elber, J. Chem. Phys. 105, 9299 (1996). 

^®R. Olender and R. Elber, J. Mol. Struct. 398, 63 (1997), 

^°C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, J. Chem. Phys. 108, 1964 
(1998). 

^^E. Weinan, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66 (2002). 

E, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden, Commun. Comput. Phys. 
(2007). 

Zhou, W. Ren, and W. E, J. Chem. Phys. 128, 104111 (2008), 

^■^A. E. Likhtman and C. M. Marques, Europhysics Letters 76, 753 (2006), 

Cao, Z. Wang, and A. E. Likhtman, to be submitted. 




